Rasters.jl

Common Applications

Subsetting an object

Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jl Dimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:

Use the .. selector to take a view of madagascar:

using Rasters, Plots
A = Raster(WorldClim{BioClim}, 5)
madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04)) # Note the space between .. -12
plot(madagascar)

Methods that change the reslolution or extent of an object

Click through to the function documentation for more in-depth descriptions and examples.

Methods that change an objects values:

Note that most regular Julia methods, such as replace, work as for a standard Array. These additional methods are commonly required in GIS applications.

Point, polygon and table operation

Methods to load, write and modify data sources:

Altering and summarising arrays and stacks with regular julia methods

Most base methods work as for regular julia Arrays, such as reverse and rotations like rotl90. Base, statistics and linear algebra methods like mean that take a dims argument can also use the dimension name. To take the mean over the time dimension:

mean(A, dims=Ti)

broadcast works lazily from disk when lazy=true, and is only applied when data is directly indexed. Adding a dot to any function will use broadcast over a Raster just like an Array.

Broadcasting

For a disk-based array A, this will only be applied when indexing occurs or when we read the array.

A .*= 2

To broadcast directly to disk, we need to open the file in write mode:

open(Raster(filename); write=true) do O
    O .*= 2
end

To broadcast over a RasterStack use map, which applies a function to the raster layers of the stack.

newstack = map(stack) do raster
    raster .* 2
end

Modifying object properties

rebuild can be used to modify the fields of an object, generating a new object (but possibly holding the same arrays or files).

If you know that a file had an incorrectly specified missing value, you can do:

rebuild(A; missingval=-9999)

(replace_missing will actualy replace the current values)

Or if you need to change the name of the layer:

rebuild(A; name=:temperature)

set can be used to modify the nested properties of an objects dimensions, that are more difficult to change with rebuild. set works on the principal that dimension properties can only be in one specific field, so we generally don't have to specify which one it is. set will also try to update anything affected by a change you make.

This will set the X axis to specify points, instead of intervals:

set(A, X => Points)

We can also reassign dimensions, here X becomes Z:

set(A, X => Z)

setcrs(A, crs) and setmappedcrs(A, crs) will set the crs value/s of an object to any GeoFormat from GeoFormatTypes.jl.

Examples and Plotting

Plots.jl is fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.

Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):

using Rasters, Plots
A = Raster(WorldClim{BioClim}, 5)
plot(A; max_res=3000)

Loading and plotting data

Our first example simply loads a file from disk and plots it.

This netcdf file only has one layer, if it has more we could use RasterStack instead.

using Rasters, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)

180×170×24 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{CFTime.DateTime360Day} CFTime.DateTime360Day[CFTime.DateTime360Day(2001-01-16T00:00:00), …, CFTime.DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (CFTime.DateTime360Day(2001-01-01T00:00:00), CFTime.DateTime360Day(2003-01-01T00:00:00)))missingval: missing
crs: EPSG:4326
mappedcrs: EPSG:4326
values: [:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.43   271.437  271.445  271.459
   3.0     missing     missing     271.431  271.438  271.445  271.459
   5.0     missing     missing     271.431  271.438  271.445  271.459
   7.0     missing     missing     271.431  271.439  271.446  271.459
   9.0     missing     missing  …  271.431  271.439  271.446  271.459
   ⋮                            ⋱                      ⋮      
 349.0     missing     missing     271.43   271.435  271.445  271.459
 351.0     missing     missing     271.43   271.435  271.445  271.459
 353.0     missing     missing     271.43   271.436  271.445  271.459
 355.0     missing     missing  …  271.43   271.436  271.445  271.459
 357.0     missing     missing     271.43   271.437  271.445  271.459
 359.0     missing     missing     271.43   271.437  271.445  271.459
[and 23 more slices...]

Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:

A[Ti=1:3:12] |> plot

Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.

using CFTime
A[Ti(Near(DateTime360Day(2001, 01, 17))), Y(-60.0 .. 90.0), X(45.0 .. 190.0)] |> plot

Now get the mean over the timespan, then save it to disk, and plot it as a filled contour:

Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.

using Statistics
# Take the mean
mean_tos = mean(A; dims=Ti)

180×170×1 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{CFTime.DateTime360Day} CFTime.DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):CFTime.DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (CFTime.DateTime360Day(2001-01-01T00:00:00), CFTime.DateTime360Day(2003-01-01T00:00:00)))missingval: missing
crs: EPSG:4326
mappedcrs: EPSG:4326
values: [:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.427  271.434  271.443  271.454
   3.0     missing     missing     271.427  271.434  271.443  271.454
   5.0     missing     missing     271.427  271.434  271.443  271.454
   7.0     missing     missing     271.427  271.435  271.444  271.454
   9.0     missing     missing  …  271.427  271.435  271.444  271.454
   ⋮                            ⋱                      ⋮      
 349.0     missing     missing     271.427  271.433  271.443  271.454
 351.0     missing     missing     271.427  271.433  271.443  271.454
 353.0     missing     missing     271.427  271.433  271.443  271.454
 355.0     missing     missing  …  271.427  271.433  271.443  271.454
 357.0     missing     missing     271.427  271.433  271.443  271.454
 359.0     missing     missing     271.427  271.433  271.443  271.454

Plot a contour plot

contourf(mean_tos; dpi=300, size=(800, 400))

Write the mean values to disk

write("mean_tos.nc", mean_tos)

"mean_tos.nc"

Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :

A[Y(Near(20.0)), Ti(1)] |> plot

A basic species distribution modelling workflow

Load occurrences for the Mountain Pygmy Possum using GBIF.jl

using Rasters, GBIF2, Plots
records = GBIF2.occurrence_search("Burramys parvus"; limit=300)

300-element GBIF2.Table{GBIF2.Occurrence, JSON3.Array{JSON3.Object, Vector{UInt8}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}}}┌──────────────────┬─────────────────┬─────────┬─────────┬─────────┬────────────
│ decimalLongitude │ decimalLatitude │    year │   month │     day │  kingdom  ⋯
│         Float64? │        Float64? │  Int64? │  Int64? │  Int64? │  String?  ⋯
├──────────────────┼─────────────────┼─────────┼─────────┼─────────┼────────────
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.327 │        -36.2604 │    2016 │      11 │      15 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│           148.22 │        -36.3511 │    2015 │      11 │      15 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.215 │        -36.4215 │    2011 │      11 │      21 │ Animalia  ⋯
│          missing │         missing │    2021 │       1 │       6 │ Animalia  ⋯
│          147.096 │        -36.9357 │    2020 │       2 │      10 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.347 │        -36.5047 │    2012 │      11 │      22 │ Animalia  ⋯
│          148.395 │        -36.5505 │    2016 │       1 │       3 │ Animalia  ⋯
│          148.236 │        -36.5249 │    2012 │      11 │      23 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│        ⋮         │        ⋮        │    ⋮    │    ⋮    │    ⋮    │    ⋮      ⋱
└──────────────────┴─────────────────┴─────────┴─────────┴─────────┴────────────
                                                 77 columns and 285 rows omitted

Extract the longitude/latitude value to a Vector of points (a Tuple counts as a (x, y) point in GeoInterface.jl):

coords = [(r.decimalLongitude, r.decimalLatitude) for r in records if !ismissing(r.decimalLatitude)]

259-element Vector{Tuple{Float64, Float64}}:
 (148.326515, -36.260389)
 (148.2198, -36.351113)
 (148.214924, -36.421454)
 (147.096394, -36.935687)
 (148.347186, -36.504673)
 (148.395217, -36.550512)
 (148.235596, -36.524924)
 (148.240881, -36.400058)
 (151.250866, -33.831883)
 (147.1, -37.0)
 ⋮
 (148.378305, -36.413602)
 (148.391682, -36.373215)
 (148.368407, -36.413914)
 (148.307214, -36.452732)
 (148.256519, -36.442973)
 (148.461854, -36.009001)
 (148.461854, -36.009001)
 (148.461854, -36.009001)
 (148.40889, -36.006805)

Get BioClim layers and subset to south-east Australia

A = RasterStack(WorldClim{BioClim}, (1, 3, 7, 12))
se_aus = A[X(138 .. 155), Y(-40 .. -25), Band(1)]

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(138.0, 154.833, 102) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(-25.1667, -39.8333, 89) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :bio1  Float32 dims: X, Y (102×89)
  :bio3  Float32 dims: X, Y (102×89)
  :bio7  Float32 dims: X, Y (102×89)
  :bio12 Float32 dims: X, Y (102×89)

Plot BioClim predictors and scatter occurrence points on all subplots

p = plot(se_aus);
kw = (legend=:none, opacity=0.5, markershape=:cross, markercolor=:black)
foreach(i -> scatter!(p, coords; subplot=i, kw...), 1:4)
display(p)

Then extract predictor variables and write to CSV.

using CSV
predictors = collect(extract(se_aus, coords))
CSV.write("burramys_parvus_predictors.csv", predictors)

"burramys_parvus_predictors.csv"

Or convert them to a DataFrame.

using DataFrames
df = DataFrame(predictors)
df[1:5, :]

Polygon masking, mosaic and plot

In this example we will mask the Scandinavian countries with border polygons, then mosaic together to make a single plot.

First, get the country boundary shape files using GADM.jl.

using Rasters, Shapefile, Plots, Dates, Downloads

# Download the shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
Downloads.download(shapefile_url, shapefile_name)

# Load using Shapefile.jl
shapes = Shapefile.Handle(shapefile_name)
denmark_border = shapes.shapes[71]
norway_border = shapes.shapes[53]
sweden_border = shapes.shapes[54]

Polygon(4665 Points)

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl, and drop the Band dimension.

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)[Band(1)]

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (2160×1080)
  :tmax Float32 dims: X, Y (2160×1080)
  :prec Int16 dims: X, Y (2160×1080)
  :wind Float32 dims: X, Y (2160×1080)

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (98×102)
  :tmax Float32 dims: X, Y (98×102)
  :prec Int16 dims: X, Y (98×102)
  :wind Float32 dims: X, Y (98×102)

Plotting

First define a function to add borders to all subplots.

function borders!(p, poly)
    for i in 1:length(p)
        plot!(p, poly; subplot=i, fillalpha=0, linewidth=0.6)
    end
    return p
end

borders! (generic function with 1 method)

Now we can plot the individual countries.

dp = plot(denmark)
borders!(dp, denmark_border)

sp = plot(sweden)
borders!(sp, sweden_border)

np = plot(norway)
borders!(np, norway_border)

The Norway shape includes a lot of islands. Lets crop them out using .. intervals:

norway_region = climate[X(0..40), Y(55..73)]
plot(norway_region)

And mask it with the border again:

norway = mask_trim(norway_region, norway_border)
np = plot(norway)
borders!(np, norway_border)

Now we can combine the countries into a single raster using mosaic. first will take the first value if/when there is an overlap.

scandinavia = mosaic(first, denmark, norway, sweden)

RasterStack with dimensions: 
  X Projected{Float64} 3.1666666666666443:0.16666666666666666:32.49999999999998 ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} 72.66666666666666:-0.16666666666666666:52.99999999999999 ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (177×119)
  :tmax Float32 dims: X, Y (177×119)
  :prec Int16 dims: X, Y (177×119)
  :wind Float32 dims: X, Y (177×119)

And plot scandinavia, with all borders included:

p = plot(scandinavia)
borders!(p, denmark_border)
borders!(p, norway_border)
borders!(p, sweden_border)

And save to netcdf - a single multi-layered file, and tif, which will write a file for each stack layer.

write("scandinavia.nc", scandinavia)
write("scandinavia.tif", scandinavia)

("scandinavia_tmin.tif", "scandinavia_tmax.tif", "scandinavia_prec.tif", "scandinavia_wind.tif")

Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.

Objects

Raster

Spatial raster data is essentially just an Array. But Raster wrappers allow treating them as an array that maintains its spatial index, crs and other metadata through all transformations. This means they can always be plotted and written to disk after applying most base Julia methods, and most broadcasts.

AbstractRaster <: DimensionalData.AbstractDimArray

Raster <: AbsractRaster

Raster(filepath::AbstractString, dims; kw...)
Raster(A::AbstractArray{T,N}, dims; kw...)
Raster(A::AbstractRaster; kw...)

Raster(T::Type{<:RasterDataSource}, [layer]; kw...) => Raster

RasterStack

Spatial data often comes as a bundle of multiple named arrays, as in netcdf. RasterStack can represent this, or multiple files organised in a similar way.

AbstractRasterStack

RasterStack <: AbstrackRasterStack

RasterStack(data...; name, kw...)
RasterStack(data::Union{Vector,Tuple}; name, kw...)
RasterStack(data::NamedTuple; kw...))
RasterStack(s::AbstractRasterStack; kw...)
RasterStack(s::AbstractRaster; layersfrom=Band, kw...)
RasterStack(filename::AbstractString; kw...)

files = (temp="temp.tif", pressure="pressure.tif", relhum="relhum.tif")
stack = RasterStack(files; mappedcrs=EPSG(4326))
stack[:relhum][Lat(Contains(-37), Lon(Contains(144))

RasterStack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => RasterStack

RasterSeries

A series is a meta-array that holds other files/data that is distributed over some dimension, often time. These files/data can be Rasters or RasterStacks.

AbstractRasterSeries <: DimensionalData.AbstractDimensionalArray

RasterSeries[Time(Near(DateTime(2001, 1))][:temp][Y(Between(70, 150)), X(Between(-20,20))] |> plot`

RasterSeries <: AbstractRasterSeries

RasterSeries(arrays::AbstractArray{<:AbstractRaster}, dims; kw...)
RasterSeries(stacks::AbstractArray{<:AbstractRasterStack}, dims; kw...)
RasterSeries(filepaths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
RasterSeries(dirpath:::AbstractString, dims; ext, child, duplicate_first, kw...)

RasterSeries(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractRasterSeries

Dimensions

Rasters uses X, Y, and Z dimensions from DimensionalData.jl to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData.jl for more details on how they work.

Band <: Dimension

Band(val=:)

banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)

Lookup Arrays

These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.

AbstractProjected <: AbstractSampled

Projected <: AbstractProjected

Projected(order, span, sampling, crs, mappedcrs)
Projected(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs, mappedcrs=nothing)

GDALarray(filename; mappedcrs=EPSG(4326))

Mapped <: AbstractProjected

Mapped(order, span, sampling, crs, mappedcrs)
Mapped(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs=nothing, mappedcrs)

Data sources

Rasters.jl uses a number of backends to load raster data. Raster, RasterStack and RasterSeries will detect which backend to use for you, automatically.

GRD

R GRD files can be loaded natively, using Julias MMap - which means they are very fast, but are not compressed. They are always 3 dimensional, and have Y, X and Band dimensions.

NetCDF

NetCDF .nc files are loaded using NCDatasets.jl. Layers from files can be loaded as Raster("filename.nc"; key=:layername). Without key the first layer is used. RasterStack("filename.nc") will use all netcdf variables in the file that are not dimensions as layers.

NetCDF layers can have arbitrary dimensions. Known, common dimension names are converted to X, Y Z, and Ti, otherwise Dim{:layername} is used. Layers in the same file may also have different dimensions.

NetCDF files still have issues loading directly from disk for some operations. Using read(ncstack) may help.

GDAL

All files GDAL can access, such as .tiff and .asc files, can be loaded, using ArchGDAL.jl. These are generally best loaded as Raster("filename.tif"), but can be loaded as RasterStack("filename.tif"; layersfrom=Band), taking layers from the Band dimension, which is also the default.

SMAP

The Soil Moisture Active-Passive dataset provides global layers of soil moisture, temperature and other related data, in a custom HDF5 format. Layers are always 2 dimensional, with Y and X dimensions.

These can be loaded as multi-layered RasterStack("filename.h5"). Individual layers can be loaded as Raster("filename.h5"; key=:layerkey), without key the first layer is used.

smapseries(filenames::AbstractString; kw...)
smapseries(filenames::Vector{<:AbstractString}, dims=nothing; kw...)

Writing file formats to disk

Files can be written to disk in all formats other than SMAP HDF5 using write("filename.ext", A). See the docs for write. They can (with some caveats) be written to different formats than they were loaded in as, providing file-type conversion for spatial data.

Some metadata may be lost in formats that store little metadata, or where metadata conversion has not been completely implemented.

RasterDataSources.jl integration

RasterDataSources.jl standardises the download of common raster data sources, with a focus on datasets used in ecology and the environmental sciences. RasterDataSources.jl is tightly integrated into Rasters.jl, so that datsets and keywords can be used directly to download and load data as a Raster, RasterStack, or RasterSeries.

using Rasters, Plots, Dates
A = Raster(WorldClim{Climate}, :tavg; month=June)
plot(A)

See the docs for Raster, RasterStack and RasterSeries, and the docs for RasterDataSources.getraster for syntax to specify various data sources.

Exported functions

Rasters.jl is a direct extension of DimensionalData.jl. See DimensionalData.jl docs for the majority of types and functions that can be used in Rasters.jl.

Functions more specific to geospatial data are included in Rasters.jl, and listed below.

aggregate(method, object, scale; filename, progress, skipmissing)

using Rasters, Statistics, Plots
using Rasters: Center
st = read(RasterStack(WorldClim{Climate}; month=1))
ag = aggregate(Center(), st, (Y(20), X(20)); skipmissingval=true, progress=false)
plot(ag)
savefig("build/aggregate_example.png"); nothing
# output